clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================


filetosave='Vilcoal_30_4_18_vill_cd_incproc_outl'; % name of workspace saved

%====================================================
% select parameters
% ===================================================

% ======== a. individual's income process ========

village=2;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=0;                         % set to 1 for coalitional deviations
outlier=0;                         % only applies to village 1: set this to 0 to load datasheet without outlier correction, 1 with larges observations eliminated, 2 with first period eliminated altogether
ytaxrate=0;                         % linear tax on income

Tessa1Tobi2=1;
n_rhogrid=25;
Nphi=10;
server=0;				% set to 1 to run on server; o wise path for workstations is used
buildup=1;              %%%set to 1 and change Qrho and Qbeta loop to calculate for all group sizes
if village==1
    vss_high_inddev=34-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[3 vss_high_inddev]; % if you choose buildup==1, rhis is the vector of village sizes for which the solution is calculated in the individual deviations case
if server==1
outputpath=sprintf('/home/workgroups/testob/Output estimation 5 December full grid outlier %d',outlier);
programpath=sprintf('/home/workgroups/testob');
else
    if Tessa1Tobi2==2
outputpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/JEEA Replication/Programmes')
programpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessaJEEA Replication/Programmes')
    elseif Tessa1Tobi2==1
        
  outputpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
programpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
      
    end
end

      betavec0=[0.9 0.94];
      betavec1=[0.95:0.01:0.96];
      betavec3=[0.94:0.01:0.98];
      if coaldev<2 || coaldev==3
      rhovec=[1];
      elseif coaldev==2
         rhovec=ones(1,Nphi); 
      end

%main execution options
gapcritset=0.01; %convergence criterion for change in values

% other execution options
Nincproc=1;                         % number of income processes to look at
T=1200; Tinit=201;                 % number of simulation periods in total, and number of initial periods to discard
maxcerr=0.9;                          % max cons meas error in log levels
N_cerr=1; %%%remember to set this back                           % number of support points for cons meas error
maxincerr=2/3;                      % max inc meas error in multiples of inc variance (in levels)
% set to one if you want to specify the rest of village in terms of average
% consumption
rov_average=1;
% set this to 1 if you want to calculate the moments for levels of
% consumption based on log consumption
momentclevinlog=1;
Rov_sb=2; % set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
% just the utility from consuming average ROV income (first-best)
% convergence criterion and maximum number of iterations
maxitgapset=400;
% this is a punishment factor: default involves a consumption penalty in the first period equal
% to caut*Pun_fac
Pun_fac=0.0;
% grid of time-varying planner weights
%lambda=[0.1,0.125,0.15,0.175,0.2,0.25,0.3,0.35:0.05:0.65,0.7,0.75,0.8,0.825,0.85,0.875,0.9]';
lambda=[0.05:0.001:0.95]';%linspace(0.05,0.95,301)';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);                  % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% I. Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, secdon: Kanzara, third: Shirapur
if outlier==0
    datafile='armoments_fe0';
elseif outlier==1
    datafile='armoments1';
elseif outlier==2
    datafile='armoments2';
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);


varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
rho_max=covyylag(village)/(vary(village)-varnu_max(village));          % AR(1) coefficient w max meas error
rho1vec=linspace(rho_min,rho_max,Nincproc);
varnu=max(vary(village)-covyylag(village)./rho1vec,0); % max prevents the first entry to be -e17
vareps=(vary(village)-varnu).*(1-rho1vec.^2);

% ===================================
% iterate over income processes
% ===================================
disp('now solving for village,coaldev outlier - press any key to continue')
disp([village coaldev outlier])
pause

%Momentsest=nan(28,11,size(betavec,2),size(rhovec,2),vss_high_inddev);
indicator1=zeros(1,size(rhovec,2));
incproc=1;   
filetosaveest=['Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(ytaxrate*10)];
% specify a markov process for individual income
[incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);


incomevals=exp(incomevals');
QQQ=Q1^1000;
SS=QQQ(1,:)';   
incomevals=incomevals/(SS'*incomevals);
cumQ=cumsum(Q1,2);
incomevals=(1-ytaxrate)*incomevals+ytaxrate; % tax as in Krueger and Perri JET 2011


% set minimum and maximum village size
if coaldev==1                       % for caolitional deviations,
   
    vss_set=[1:12,17,22,27,vss_high_coaldev]
    vss_high=vss_set(end);
elseif coaldev==0
    if buildup==0% for individual deviations, start with small village and build up, but coarser
        vss_set=[vss_high_inddev];
        lagind=1;
        vss_high=vss_set(end);% an indicator
    elseif buildup==1
        vss_set=vssbuildupvec_inddev;
        lagind=1;
        vss_high=vss_set(end);% an indicator
    end
elseif coaldev==2;
    vss_set=vss_high_inddev;
     vss_high=vss_set(end);
elseif coaldev==3
    vss_set=[1:10];
     vss_high=10;
end


% ===================================
% iterate over discount factors
% ===================================
if coaldev==0 
            betavec=betavec0;
elseif coaldev==1
               betavec=betavec1; 
elseif coaldev==2
    r=0.04;
    betamax=1/(1+r)-0.015;
    betavec=linspace(0.8,betamax,16);
             Z=incomevals;
          phinat=min(Z)/(r-1);
          phivec=linspace(phinat,0,Nphi);
elseif  coaldev==3
    betavec=betavec3; 
end
            Income0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nincproc);
Csim0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nincproc);
Momentsest=nan(30,N_cerr+1,size(betavec,2),1,vss_high_inddev);
for Qbeta=1:size(betavec,2) 
    
    % ===================================
    % iterate over CRRA
    % ===================================
    for Qrho=1:size(rhovec,2) 
beta=betavec(Qbeta)
rho=rhovec(Qrho)  
            
            % limited commitment economies
            
        if coaldev==1
            lagind=[];
        else
            lagind=1;
        end
       
        P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
        waut_lag=[]; waut_avg_lag=[];
         WAUT_LAG=nan(1000,1000,length(vss_set));
            STABLE=[];
        for vsscount=1:length(vss_set)       % iterate over village size
            vss=vss_set(vsscount);
            vss_max=max(vss_set);
            if vsscount>1 & coaldev==1
                lagind=lagind+(vss-(vss_set(vsscount-1)+1));
            end
            
            disp('now solving for parameters')
            disp([Qbeta,Qrho])
            disp('vss')
            disp(vss)
             % Aiyagari economy     
        if coaldev==2
            sav_problem
            % Full insurance with coal deviations - simple example in
            % Genicot and Ray
        
        else
            % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
            % by vss
            if rov_average==1
                ScaleRov=vss;
            else
                ScaleRov=1;
            end
            csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; waut_check=[]; waut_check_ind=[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
            S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
            wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
            Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
            S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all=[];S_rov_aut_all2=[];
            w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];
            
            
            % ===================================
            % build the state space of ROV income
            % ===================================
            idio_dist=[];
            idio_dist_aut=[];
            idio_dist_aut_all=[];idio_dist_aut_all2=[];
            idio_dist_sb=[];
            
            for k1=0:vss
                for k2=0:vss
                    for k3=0:vss
                        II=zeros(1,3);
                        if k1+k2+k3==vss
                            % the vector of aggregate incomes from
                            % permutations of vss households
                            S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1) + k3*incomevals(3,1)];
                            % the corresponding distribution of
                            % individual income values
                            idio_dist=[idio_dist; k1 k2 k3];
                            
                            % this computes the 'best' distribution of incomes
                            % of size vss-lagind
                            if vss>1
                                II=zeros(1,3);
                                temp=[k1 k2 k3];
                                % this takes away the lowest incomes to
                                % derive the distribution of incomes of the
                                % 'best' coalition lagind individuals
                                % This coalition is one smaller than the
                                % largest sustainable coalition because it
                                % will consist of 1 agent and m-1 guys from
                                % the largest sustainable coalition of size
                                % m
                                
                                for jj=1:lagind
                                    II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                    temp=[k1 k2 k3]-II;
                                end
                                idio_dist_aut=[idio_dist_aut; [k1 k2 k3]-II];
                                S_rov_aut=[S_rov_aut; ([k1 k2 k3]-II)*incomevals];
                                
                                                               
                                II=zeros(1,3);
                                temp=[k1 k2 k3];
                                % this takes away the lowest incomes to
                                % derive the distribution of incomes of the
                                % 'best' coalition
                                % consisting of m individuals from only the
                                % rest of village
                                for jj=1:lagind-1
                                    II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                    temp=[k1 k2 k3]-II;
                                end
                                idio_dist_sb=[idio_dist_sb; [k1 k2 k3]-II];
                                S_rov_sb=[S_rov_sb; ([k1 k2 k3]-II)*incomevals];
                            end
                        end
                    end
                end
            end
            
            % resort S_rov in ascending order
            [S_rov,S_rovind]=sort(S_rov);
            idio_dist=idio_dist(S_rovind,:);
            
            if vss>1
                idio_dist_aut=idio_dist_aut(S_rovind,:);
                S_rov_aut=S_rov_aut(S_rovind);
                idio_dist_sb=idio_dist_sb(S_rovind,:);
                S_rov_sb=S_rov_sb(S_rovind);
            end
            
            
            if vss>1
            for k=1:length(STABLE)
            if STABLE(k)==1 
            for g=1:length(idio_dist)
            index=0;
            for j1=0:vss-k 
                for j2=0:vss-k
                    for j3=0:vss-k
                        if j1+j2+j3==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3)
                           index=index+1  ;   
                           temp=[j1 j2 j3];
                                idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp]; 
                                S_rov_aut_all(g,:,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                        end      
                     end
                 end
            end
            end                    
            end
            
            end
            end
            % ===================================
            % simulate aggregate income transitions
            % ===================================
            tic
            T_trans=50000;
            aggincind=zeros(length(S_rov),length(S_rov));
            for jjj=1:length(S_rov)
                jjj;
                inccount=zeros(T_trans,vss,3);
                agginc=zeros(T_trans,3);
                for iii=find(idio_dist(jjj,:)>0)
                    
                    for kkk=1:idio_dist(jjj,iii)
                        rng(jjj*100+iii*10+kkk,'twister');
                        shock=rand(T_trans,1);
                        rr=(find(shock<=cumQ(iii,1)));
                        inccount(rr,kkk,iii)=1;
                        rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                        inccount(rr,kkk,iii)=2;
                        rr=(find(cumQ(iii,2)<shock & cumQ(iii,3)>=shock));
                        inccount(rr,kkk,iii)=3;
                    end
                    
                    agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                end
                agginc=sum(agginc,2);
                for jjj2=1:length(S_rov)
                    qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                    aggincind(jjj,jjj2)=sum(qqq);
                end
                P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
            end
            disp('Simulation for agg income process takes')
            toc
            
            
           if vss>1
                S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
                S_rov_sb=S_rov_sb/(max(ScaleRov-lagind+1,1));
                for g=vss_set(find(vss-1>=vss_set))
                if STABLE(g)==1
                S_rov_aut_all(:,:,:,g)=S_rov_aut_all(:,:,:,g)/(g);
                end
                end
            end
            
            S_rov=S_rov/ScaleRov;
            
            % ===================================
            % state space: combinations of indiv and village income
            % ===================================
            S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
            idio_dist_rov=kron(ones(size(incomevals)),idio_dist);
            idio_dist_aut=kron(ones(size(incomevals)),idio_dist_aut);
            S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
            S_rov_sb=[kron(ones(size(incomevals)),S_rov_sb)];
            
            if vss>1
            for g=vss_set(find(vss-1>=vss_set))
            if STABLE(g)==1    
            for k=1:size(S_rov_aut_all,3)
            S_rov_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),S_rov_aut_all(:,:,k,g))];
            idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),idio_dist_aut_all(:,:,k,g))];
            end
            end
            end
            S_rov_aut_all=S_rov_aut_all2;
            idio_dist_aut_all=idio_dist_aut_all2;
            clear S_rov_aut_all2;
            clear idio_dist_aut_all2;
            end
            
            P=kron(Q1,P_rov);
            
            Y=S(:,1)+ScaleRov*S(:,2);				% possible aggregate income outcome
            state=length(S);                    % number of states
            agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
            STATE(vss)=state;
            
            
            % ===================================
            % Compute set of autarky values
            % ===================================
            for k=1:agent
                caut(:,k)=S(:,k);
            end
            waut=zeros(state,agent);
            waut_avg_lag=zeros(state,agent);
            waut_sb=zeros(state,agent);
            % this is just the PDV of utility from consumption of
            % income caut
            eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
            eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
            if coaldev<2
            % in first iteration, or for individual deviations, waut is just consumption of income
            if vss==1 || coaldev==0
                for j=1:state
                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                    % autarky values for the rest of village are always the same - perfect
                    % risk-sharing among the rest of village
                    waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                end
                
            else
                % this calculates autarky values for the coalitional
                % deviations case
                               
                for j=1:state
                    
                    S_lag_ind_sb=find((S_rov_sb(j,1)<=(S_lag(:,1)+S_lag(:,2)*(vss-lagind))/(vss-lagind+1)+0.00001 & S_rov_sb(j,1)>=(S_lag(:,1)+S_lag(:,2)*(vss-lagind))/(vss-lagind+1)-0.00001));
                   
                    wwww=[utility(rho,S_lag(S_lag_ind_sb,1))+beta*P_lag(S_lag_ind_sb,:)*waut_lag(:,1) utility(rho,S_lag(S_lag_ind_sb,2))+beta*P_lag(S_lag_ind_sb,:)*waut_lag(:,2)];
                    
                    wwwalt(j)=mean(mean(wwww,2));
                    www=(waut_lag(S_lag_ind_sb,1)+(vss-lagind)*waut_lag(S_lag_ind_sb,2))/(vss-lagind+1);
                    
                    SSS =P_lag^1000;
                    SSSS=SSS(:,1); % ergodic prob distribution
                    waut_avg_lag(j,2)=SSSS(S_lag_ind_sb)'/sum(SSSS(S_lag_ind_sb))*www;
                end
                
                
                for j=1:state
                    % indicator for the 'autarky' state in the state
                    % space of the previous iteration S_lag
                    S_lag_2=(S_rov_aut(j)==S_lag(:,2)) ;
                    S_lag_1=(S(j,1)==S_lag(:,1));
                    S_lag_ind=find(S_lag_1.*S_lag_2==1);
                    % for the individual, autarky value is utility of
                    % current consumption plus relevant trans
                    % probability times waut_lag, the vector of values
                    % from the next smaller insurance group
                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                    % for rest of village, the autarky value is just
                    % expected disc utility of current consumption
                    waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                    
                    for g=1:size(S_rov_aut_all,4)
                        if STABLE(g)==1
                        for k=1:size(S_rov_aut_all,3)
                    S_lag_2=(S_rov_aut_all(j,1,k,g)==S_LAG(:,2,g)) ;
                    S_lag_1=(S(j,1)==S_LAG(:,1,g));
                    S_lag_ind=find(S_lag_1.*S_lag_2==1);
                        
                    agentstate=[];
                    for kk=1:length(incomevals)
                        if idio_dist_aut_all(j,kk,k,g)>0
                    agentstate=[agentstate kk];    
                        end
                    end
                    
                    if size(S_lag_ind)>0
                    for kk=1:length(agentstate) %%%note, you only need one entry per income state, not one entry per agent!!!!!
                    waut_check_ind(j,1,k,g)=utility(rho,S(j,1))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),1,g);
                    waut_check(j,kk,k,g)=utility(rho,incomevals(agentstate(kk)))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),2,g);
                    S_check(j,kk,k,g)=agentstate(kk);  %%%not sure that we need this
                    
                    end
                    end
                    
                        end
                    end
                    end
                    
                    
                    
                    %%%and the second best utilities
                    waut_sb(j,1)=waut(j,1);
                    waut_sb(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_avg_lag(:,2);
                    waut_sb2(j,1)=waut(j,1);
                    waut_sb2(j,2)=wwwalt(j);
                    
                if S(j,1)==incomevals(3) && idio_dist_rov(state,3)==vss-1 
                   waut_alt(j,1)=waut(j,:);
                else 
                    waut_alt(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                   
                    waut_alt(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                end
                end

            end
            
                     
            
            if Rov_sb==1 & vss>1 & coaldev==1
                waut=waut_sb;
            end
            
            % ====================================================
            % Compute constrained insurance
            % ====================================================
            gammagrid=[lambda vss/ScaleRov*(1-lambda)];
            gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
	        gammagridold=gammagrid;
	        clear gammagrid
	        gammagrid=[linspace(1.2*max(S(:,1)),min(S(:,1)),m)',linspace(min(S(:,2)),1.2*max(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,1)),min(S(:,1)),m)'./linspace(min(S(:,2)),1.2*max(S(:,2)),m)'];
            
            % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
            n=length(gammagrid);
            w=zeros(n,state,agent);
            gammar=zeros(n,state,agent);
            c=zeros(n,state,agent);
            phi=zeros(n,state,agent);
            
            % 'first-best' consumption given aggregate resources and planner weights
            for j=1:state
                c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
            end
            
            cfirstbest=c;
            %  relative planner weights
            for j=1:state
                for k=1:agent
                    gammar(:,j,k)=gammagrid(:,agent+k);
                end
            end
            
            % Compute the first best value functions without binding
            % constraints - i.e. cfirstbest for ever
            w=zeros(n,state,agent);
            
            dww=1;
            wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
            while dww>0.00001
                w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                
                dww=max(max(max(abs(w-wold))));
                wold=w;
            end
            wfirstbest=w;
            
            % continuation utility as function of planner weight and state of the
            % economy
            ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
            ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';
            
            
            % ====================================================
            % Start the loop with the enforcement constraints
            % ====================================================n
            ewsb=ewfirstbest;
            wsb=wfirstbest;
            options=optimset('Display','off');
            eps=0.0000001;t=0;itgap=0; gap=2; gaphist=gap; gapincrease=0;
            
            csolold=cfirstbest;
	     if beta<0.85
		gapcrit=gapcritset;
		maxitgap=maxitgapset;
		else
		gapcrit=gapcritset*1
		maxitgap=maxitgapset;
         end
        
        
            while gap >=gapcrit && itgap<=maxitgap
                rr=[];
                rr1=[];
                rr2=[];
                rrnone=[];
                rrboth=[];
                nonconv=zeros(state,3)*nan;
                itgap=itgap+1;
                coalnonsust=nan(state,2);
                count=0; countgrid=[ones(n,state)];
                for j=1:state
                                  
                    %%instead of solving problem for each grid point,
                    %%note that there are only three possibilities: 1)
                    %%Constraint binding for Person 1 only, the
                    %%solution is the same for each lambda for which
                    %%this is the case
                    
                    index=[1 j];
                    
                    rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                    if vss>1 & coaldev==1
                       rr1cd=[]; 
                        for g=size(waut_check,4):-1:0 %%%check from the largest group down, 0 is an individual deviation
                        
                        if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                        %%%check the largest possible coalition first:
                        if g>0
                        for k=1:size(waut_check,3)  
                            if waut_check(j,1,k,g)~=0
                            coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-waut_check_ind(j,1,k,g)]; 
                            for kk=1:size(waut_check,2)
                                if waut_check(j,kk,k,g)~=0
                            coalition=[coalition utility(rho,cfirstbest(rr1,j,2))+beta*ewsb(rr1,j,2)-waut_check(j,kk,k,g)];
                                end
                            end
                       %%%Step 1: find the grid points for which you are now trying to find a solution  
                       
                       %%%%do the others want to come along?
                       rrcheck=coalition<0;    
                       rrcheck=prod(rrcheck,2);
                       rrcheck=find(rrcheck>0);
                       if size(rrcheck>0)
                       
                       bind=1;
                        
                       %%%%Step 2: find the solution
                       coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-waut_check_ind(j,1,k,g)]; 
                            for kk=1:size(waut_check,2)
                                if waut_check(j,kk,k,g)~=0
                            coalitionsolve=[coalitionsolve utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)-waut_check(j,kk,k,g)];
                                end
                            end
                       rrchecksolve=coalitionsolve<0;    
                       rrchecksolve=prod(rrchecksolve,2);
                       rrchecksolve=find(rrchecksolve>0);  
                       
                       rr=rrchecksolve(end)+1;
                            if rr>length(gammagrid)
                                coalnonsust(j,1)=-2;
                                output=[nan,nan,-50];
                        else
                            
                            clear x0
                            
                            x0(1) = cfirstbest(rr,j,1);  
                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[waut_check_ind(:,1,k,g) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                                           
                            
                       end
                            if output(3) <=0
                                nonconv(j,1)=output(3);
                                csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                
                                                               
                            else
                                
                               
                                csol(rr1(rrcheck),j,1)=x(1);
                                gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                phisol(rr1(rrcheck),j,2:agent)=0	;
                            end
                            end
                            end
                        end
                        
                        
                        rr1(rrcheck)=[];  %%%these grid points have been dealt with
                        rrnoneadd=rr1;
                        elseif g==0 %%%only an individual deviation
                        coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-eeewaut(j,1)];     
                      
                            
                       
                       rrcheck=coalition<0;    
                       rrcheck=prod(rrcheck,2);
                       rrcheck=find(rrcheck>0);
                       if size(rrcheck>0)
                          
                       bind=1;
                        coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-eeewaut(j,1)];     
                      
                       
                       rrchecksolve=coalitionsolve<0;    
                       rrchecksolve=prod(rrchecksolve,2);
                       rrchecksolve=find(rrchecksolve>0);
                       
                       
                        
                        rr=rrchecksolve(end)+1;
                        
                        if rr>length(gammagrid)
                                coalnonsust(j,1)=-2;
                                output=[nan,nan,-50];
                        else
                            
                            clear x0
                            
                            x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                            %T vectors of consumption and utilty that meet
                            %participation constraints
                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                            
                            %T 29_Mar This keeps track of nonconvergence, but
                            %allows the programme to continue - sometimes it
                            %converges in a subsequent iteration
                            
                            
                       end
                            if output(3) <=0
                                nonconv(j,1)=output(3);
                                csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                
                                                               
                            else
                                
                               
                                csol(rr1(rrcheck),j,1)=x(1);
                                gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                phisol(rr1(rrcheck),j,2:agent)=0	;
                            end
                            end                     
                        rr1(rrcheck)=[];  %%%these grid points have been dealt with
                        rrnoneadd=rr1;
                        end
                         
                        end
                        end
                       else
                        %%%%in first group size iteration, i.e. for vss=1
                        rr1cd=rr1;
                        rr1id=[];
                        rrnoneadd=[];
                    end
                    
                           
                    
                    if size(rr1cd)>0
                        bind=1;
                        
                        
                        rr=rr1cd(end)+1;
                        if isempty(rr)==1
                                coalnonsust(j,1)=-2;
                                output=[nan,nan,-50];
                        else
                            if size(rr1cd)>0
                            clear x0
                            
                            x0(1) = cfirstbest(rr1cd(1),j,1); 
                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                          
                            end
                            
                       end
                            if output(3) <=0
                                nonconv(j,1)=output(3);
                                csol(rr1cd,j,1)=ones(size(rr1cd,1),1)*caut(j,1);
                                csol(rr1cd,j,2)=ones(size(rr1cd,1),1)*caut(j,2);
                                gammar(rr1cd,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                wsbnewsol(rr1cd,j,1)=ones(size(rr1cd,1),1)*waut(j,1);
                                wsbnewsol(rr1cd,j,2)=ones(size(rr1cd,1),1)*waut(j,2);
                                phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                phisol(rr1cd,j,2:agent)=0	;
                                
                                                               
                            else
                                
                               
                                csol(rr1cd,j,1)=x(1);
                                gammar(rr1cd,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                csol(rr1cd,j,2)=(Y(j)-x(1))/ScaleRov;
                                wsbnewsol(rr1cd,j,1)=output(1);
                                wsbnewsol(rr1cd,j,2)=output(2);
                                phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                phisol(rr1cd,j,2:agent)=0	;
                            end
                    end
                    %%%%this becomes relevant for vss>1, these are all the left-overs that are not dealt with by joint deviations.    
                    if size(rrnoneadd)>0
                        
                        csol(rrnoneadd,j,1:agent)=cfirstbest(rrnoneadd,j,1:agent);
                        gammar(rrnoneadd,j,1:agent)=gammagrid(rrnoneadd,agent+1:2*agent);
                        phisol(rrnoneadd,j,1:agent)=0;
                        
                        wsbnewsol(rrnoneadd,j,1)=utility(rho,cfirstbest(rrnoneadd,j,1))+beta*ewsb(rrnoneadd,j,1);
                        wsbnewsol(rrnoneadd,j,2)=utility(rho,cfirstbest(rrnoneadd,j,2))+beta*ewsb(rrnoneadd,j,2);
                        
                        
                        
                    end
                    
                    % 2. Constraint binding for agent 2 (RoV), 
                    rr2=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ; 
                    
                    if size(rr2)>0
                        bind=2; %%i.e. agent 2 has the binding constraint
                        rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                        if isempty(rr)==1
                                coalnonsust(j,2)=-2;
                         output=[nan,nan,-50];
                        else
                            clear x0
                            x0(1) = cfirstbest(rr(end),j,1); %%redundant because we no longer need an initial guess
                          
                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut(:,1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                        
                                               
                        end
                            % output(3) equals 1 if constraints
                            % are satisfied, but -2 if not
                            if output(3)<=0
                                nonconv(j,2)=output(3)*100;
                                csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                phisol(rr2,j,1)=0;
                            else
                             
                            if vss>1  & coaldev==1 % only need to check this once coalitional deviations become relevant
                            for g=size(waut_check,4):-1:1 %%%check from the largest group down, 0 is an individual deviation
                        
                            if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                        %%%check the largest possible coalition first:
                        
                        for k=1:size(waut_check,3)   %
                            if waut_check(j,1,k,g)~=0
                            coalition=[output(1)-waut_check_ind(j,1,k,g)]; 
                            for kk=1:size(waut_check,2)
                                if waut_check(j,kk,k,g)~=0
                            coalition=[coalition output(2)-waut_check(j,kk,k,g)];
                                end
                            end
                       %%%Step 1: find the grid points for which you are now trying to find a solution  
                       
                       %%%%do the others want to come along?
                       rrcheck=coalition<0;    
                       rrcheck=prod(rrcheck,2);
                       rrcheck=find(rrcheck>0);
                       if size(rrcheck)>0
                          output(3)=-2; 
                       end
                       
                            end
                       end
                                                        
                        end
                            end
                            end
                                
                                if output(3)>0
                               
                                csol(rr2,j,1)=x(1);
                                gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                                wsbnewsol(rr2,j,1)=output(1);
                                wsbnewsol(rr2,j,2)=output(2);
                                phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                phisol(rr2,j,1)=0;
                                else
                                 nonconv(j,2)=output(3)*100;
                                csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                phisol(rr2,j,1)=0;    
                                    
                                    
                                
                                
                                end
                                end
                        
                    end
                    
                    % 3. Constraint all slack
                    rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                    
                    if size(rrnone)>0
                        
                        csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                        gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                        phisol(rrnone,j,1:agent)=0;
                        wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                        wsbnewsol(rrnone,j,2)=utility(rho,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                        
                        
                        
                    end
                    
                end
                
                
                gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
                gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                gap=max(gap1,gap2);
                gaphist(itgap+1)=gap;
                disp('gap,gapold')
                disp([gap,gapold])
                nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
		  Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
               
                             
               if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                    itgap
                    disp('Autarky or coalition not sustainable')
                    [nonconv,S]
                    for j=1:state
                       csol(:,j,1)=caut(j,1);
                       csol(:,j,2)=caut(j,2);
                    end
                   % pause
                   break
                end

                csolold=csol;
                wsb=wsbnewsol;
                  ewsb(:,:,1)=(P*wsb(:,:,1)')'; 
            end
            
              CRITGAP(vss)=gaphist(itgap-1);
               
              
            
            if itgap>=maxitgap
            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
            end
            
            if gapincrease==1
                 nonconvind(Qbeta,Qrho,vss,incproc)=-gap;
            end
            % =============================================================
            % calculate values of outside option for coalitional deviations
            % =============================================================
            % if insurance group of size vss was unsustainable, use
            % previous value for outside option
           
            if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>1
                % this is redundant, but shows you that if vss is
                % unsustainable, the outside option for vss+1 remains that
                % for vss
                waut_lag(:,1)=waut_lag(:,1);
                waut_lag(:,2)=waut_lag(:,2);
                P_lag=P_lag;
                S_lag=S_lag;
                % if the current insurance group is not sustainable, it
                % cannot be the outside option for the next bigger group,
                % so increase the 'lagind' indicator by one
                lagind=lagind+1;
                gap=-1;
                STABLE(vss)=0;
                STABLECHECK=[1 STABLE]; 
                % ... or if there is no previous, use just conditional exp of
                % autarky consumption as calculated above
            elseif (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)  && itgap>=3 && (vss==1 || coaldev==0)
                waut_lag(1:length(eeewaut),1)=eeewaut(:,1);
                waut_lag(1:length(eeewaut),2)=eeewaut(:,2);
                P_lag=P;
                S_lag=S;
                lagind=1;
                gap=-1;
                
                
                S_rov_KEEP(1:state/3,vss)=S_rov;
                S_KEEP(1:state,:,vss)=S;
                PHI_KEEP(:,1:state,:,vss)=phisol;
                STABLE(vss)=1;
                STABLECHECK=[1 STABLE]; 
                IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
            else 
                lagind=1;
                % when ins group of size vss is sustainable, its value
                % becomes outside option for next bigger group size
                for j=1:state
                    % outside option for next larger insurance group is
                    % value with equal weights vss ...
                    eewaut(j,1)=wsb((n+1)/2,j,1);
                    eewaut(j,2)=wsb((n+1)/2,j,2);
                    
                end
                waut_lag=eewaut;
                S_lag=S;
                P_lag=P;
                
                S_rov_KEEP(1:state/3,vss)=S_rov;
                S_KEEP(1:state,:,vss)=S;
                PHI_KEEP(:,1:state,:,vss)=phisol;
                WAUT_LAG(1:state,1:2,vss)=waut_lag;
                P_LAG(1:state,1:state,vss)=P_lag;
                S_LAG(1:state,:,vss)=S_lag;
                
                STABLE(vss)=1;
                STABLECHECK=[1 STABLE]; %%%this includes stability for the individual as well, makes the loop above easier
                IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
            end
            if itgap>=maxitgap
            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
            end
                      
            end
        end

        
        end   
              
        
        %==========================================================
        % Simulate the economy: CD only for the largest group size;
        % ID for all
        % =========================================================
        % This lists all possible combinations of individual income
        % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
        % with every individual income value
        
 
        
        if coaldev==1                       % for coalitional deviations, all the stable ones
            stable=(find(STABLE==1));
            vss_set_sim=stable;
        elseif coaldev==0                                % for individual deviations, buildup (but with holes)
            vss_set_sim=vss_set;
        end
        
        for vss=[vss_set_sim]
            if coaldev<2
            state=STATE(vss);
            if rov_average==1
                ScaleRov=vss;
            else
                ScaleRov=1;
            end
            
            S=S_KEEP(1:state,:,vss);
            S_rov=S_rov_KEEP(1:state/3,:,vss);
            idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss);
            phisol=PHI_KEEP(:,1:state,:,vss);
                        IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
            
            distribution=size(IDIO_DIST,1);
            end
            
            if village==1
            Nvill=ceil((vss_high_inddev+1)/(vss+1))*15  %%%this adjusts for sampling from the village, which is what you need in the correlation test later on
            elseif village==2
            Nvill=ceil((vss_high_inddev+1)/(vss+1))*6  %%%this adjusts for sampling from the village, which is what you need in the correlation test later on
            elseif village==3
            Nvill=ceil((vss_high_inddev+1)/(vss+1))*9  %%%this adjusts for sampling from the village, which is what you need in the correlation test later on
            end
            
            
            var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
            varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
            vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
            vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
            
            betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
            beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
            beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
            
            
            
            %=============================================
            % Start simulation
            % ============================================
            % number of insurance groups equals village size
            % vss_high+1 divided by size of ins group vss+1
            
            stream = RandStream.getGlobalStream;
            Npanel=Nvill*(vss+1);
            csim0=nan(Nvill,vss+1,T);
            income0=nan(Nvill,vss+1,T);
            income0ind=nan(Nvill,vss+1,T);
            phisim0=nan(Nvill,vss+1,T);
            Yvill=nan(1,T);
            Yagg=nan(Nvill,T);
            Ssim=nan(Nvill,vss+1,T);
            gammarfun=nan(Nvill,vss+1,T);
            phinew=nan(Nvill,vss+1,T);
            indSjj=nan(T,vss+1);
            
            disp('simulating')
            
            disp('Reminder: now simulating for village,coaldev,incproc')
            disp([village coaldev incproc])
            disp('now simulating for parameters')
            disp([Qbeta,Qrho])
            disp('vss')
            disp(vss)
            
            
            tic
            
            %initial state
            t=1;
            % initial values of income and consumption (set =
            % income)
            rng(3,'twister');
            shockinit=randi(3,Nvill,vss+1);
            income0ind(:,:,1)=shockinit;
            income0(:,:,1)=incomevals(income0ind(:,:,1));
            csim0(:,:,1)=income0(:,:,1);
            
            
            Yagg(:,t)=sum(income0(:,:,t),2);
            
            
            IDIO_DIST_sim=[];
            
            while t<T
                t=t+1;
                %%%%simulate next period income taking into account
                %%%% persistence across time
                rng(t,'twister');
                shock=rand(Nvill,vss+1);
                for k=1:vss+1
                    
                    rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                    income0ind(rr,k,t)=1;
                    rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                    income0ind(rr,k,t)=2;
                    rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                    income0ind(rr,k,t)=3;
                    income0(:,k,t)=incomevals(income0ind(:,k,t))  ;
                end
                
                Yagg(:,t)=sum(income0(:,:,t),2);
            end
            
            if coaldev<2
                
            t=1;
                 while t<T
                 t=t+1;
                
                gammarfun(:,:,t)=(Yagg(:,t-1)*ones(1,vss+1)-csim0(:,:,t-1))./csim0(:,:,t-1)/ScaleRov;
                gammarfun(:,:,t)=min(gammarfun(:,:,t),max(gammagrid(:,4))*ones(Nvill,vss+1));
                
                
                
                for ivill=1:Nvill
                    IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,:,t)==3))))==1)';
                    for k=1:vss+1
                        % choose for individual k the aggregate state
                        % that applies
                          indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t));
                        % check that aggregate and individual income0s are
                        % consistent with the states caculated above
                        if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                            stop
                        end
                        
                       
                    end
                    Ssim(ivill,:,t)=indSjj(t,:);
                end
                
                
                for k=1:vss+1
                    BLA=interp1(gammagrid(:,4),phisol(:,Ssim(:,k,t),1),gammarfun(:,k,t));
                    phinew(:,k,t)=diag(BLA);
                end
                
                gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
                                
                csim0(:,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(:,t)*ones(1,vss+1));
              end
            AA=[];
            BB=[];
            CC=[];
            for t=1:T
                AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
            end
            csim0keep=csim0;
            income0keep=income0;
            income0indkeep=income0ind;
                        csim0=AA;
            income0ind=CC;
            income0=BB;
            end
            disp('simulating takes')
            tsim(Qbeta,Qrho)=toc
            
            tic

           
            
            disp('reshaping takes')
            toc
            
            
            
            % =========================================================
            % calculate moments of consumption and income growth
            % =========================================================
            csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
            rng(3, 'twister'); cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
            rng(3, 'twister'); incerr=randn(size(csim0));   Momentsest1=nan(30,N_cerr); ind=[];
            % loop over size of cons measurement error
           for kkk=1:size(csim0,1)
                tempvar(kkk)=var(log(csim0(kkk,:)));
            end
            Tempvar=mean(tempvar);
            for ic=1:N_cerr 
                csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                logcsim=log(csim);
                logcsim0=log(csim0);
                logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
                % add inc meas error back in
                income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                logincome=log(income);
                logincome0=log(income0);
                logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                
                
                for ivill=1:Nvill
                    aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    logaggvillincome=log(aggvillincome);
                    
                    logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    
                    
                end
                aggincome=sum(income,1);
                
                
                
                % === discard first Tinit-1 periods
                              
                    ccc=logcsim(:,Tinit:T);
                    ccc_cond=logcsim_meandev(:,Tinit:T);
                    ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
               
                yyy=logincome(:,Tinit:T);
                yyy_cond=logincome_meandev(:,Tinit:T);
                yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                
                cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
                cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                
                yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                aggcc=log(sum(csim0(:,Tinit:T),1))-log(sum(csim0(:,Tinit-1:T-1),1));
                aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                
                
                vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                
                % calculate moments for each individiual
                for ind=1:size(csim,1)
                    
                    dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                    dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                    dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                    dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
                    
                     
                   
                    
                    
                    
                end
                 XXreg=reshape(aggvillyy,size(aggvillyy,1)*size(aggvillyy,2),1);
                    XXreg=[ones(length(XXreg),1) XXreg];
                    YYreg=reshape(cc,size(cc,1)*size(cc,2),1);
                    bb=regress(YYreg,XXreg);
                    res_cc_cond=YYreg-XXreg*bb;
                    res_cc_cond=reshape(res_cc_cond,size(cc,1),size(cc,2));  
                if buildup==1
                
                DC_LOG_COND(1:size(cc_cond,1),1:size(cc_cond,2),ic,vss,Qbeta)=cc_cond;
                DC_LOG(1:size(cc,1),1:size(cc,2),ic,vss,Qbeta)=cc;
                RES_COND(1:size(res_cc_cond,1),1:size(res_cc_cond,2),ic,vss,Qbeta)=res_cc_cond;
               
               
                end

            
            end
%             
            
            toc
            
        end
        
        % pause
      end
end

cd(outputpath)
if buildup==1 & village==1
DC_LOG_COND2=DC_LOG_COND(:,:,1,:,:);
DC_LOG_COND2=squeeze(DC_LOG_COND2);
save('permmodel1.mat','DC_LOG_COND2')
elseif buildup==1 & village==2
DC_LOG_COND2=DC_LOG_COND(:,:,1,:,:);
DC_LOG_COND2=squeeze(DC_LOG_COND2);
DC_LOG_2=DC_LOG(:,:,1,:,:);
DC_LOG_2=squeeze(DC_LOG_2);
RES_COND2=RES_COND(:,:,1,:,:);
RES_COND2=squeeze(RES_COND2);
save('permmodel2.mat','DC_LOG_COND2','DC_LOG_2','RES_COND2') 
elseif buildup==1 & village==3
DC_LOG_COND2=DC_LOG_COND(:,:,1,:,:);
DC_LOG_COND2=squeeze(DC_LOG_COND2);
save('permmodel3.mat','DC_LOG_COND2') 
end
cd(programpath)
